Set up directory, install packages and load packages

Exercise 1: Create study area map using “ggplot”

All datasets are found in the “data”-folder.

Here, we use the package ‘rnaturalearth’ to download countries in a ESRI shapefile format. These are available for download here: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-0-countries/

Coordinate reference systems (CRS)
CRS provide a standardized way of describing locations to describe geographic data. The CRS that is chosen depends on when the data was collected, the purpose of the data, etc. It is necessary to transform vector and raster data to a common CRS when data with different CRS are combined.

Setting CRS to SWEREF99 (Sweden). CRS can be referenced by its:

  1. EPSG code, CRS(“+init=epsg:3006”) (see http://www.epsg-registry.org/ and http://spatialreference.org/) or by

  2. proj4spring; “+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs”

There are two general options:

  1. unprojected (a.k.a. Geographic): Latitude/Longitude for referencing location on the ellipsoid Earth (wgs84), and

  2. projected: Easting/Northing for referencing location on 2D representations of Earth (the creation of maps) e.g., SWEREF99.

More reading:
https://www.nceas.ucsb.edu/sites/default/files/2020-04/OverviewCoordinateReferenceSystems.pdf

Load vector data

How to download using “rnaturalearth” and save the shp locally

I will open the saved .shp for the exercise:

## load country boarders -----------


#' downladed and saved from rnaturalearth (above)
swe_no_sweref <-  input_vector %>% 

  # value = T: return vector containing the matching elements
  grep(pattern = "swe_no", value = T) %>% 

  # read ESRI Shapefile object 
  sf::st_read() %>% 

  #convert geometry object into an sf object
  st_as_sf() 
## Reading layer `swe_no' from data source 
##   `/Users/hirn0001/Library/CloudStorage/OneDrive-Sverigeslantbruksuniversitet/Undervisning/Teaching and courses/GS_VMAS Data handling and illustrations/2025/create-maps/input_data/vector/country/swe_no.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -362704.5 ymin: -6097519 xmax: 1114145 ymax: 8972413
## Projected CRS: SWEREF99 TM
# check coordinate system 
st_crs(swe_no_sweref)$proj4string
## [1] "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
# Create a bounding box as an sf object for layers ----
b_box <- c(xmin = 350432, ymin = 7179284, 
           xmax = 857780, ymax =7444230)

# Now intersect the full world's coastline with the bounding box 
suppressWarnings({ 
sweden.c <- st_crop(swe_no_sweref, b_box)
plot(sweden.c)
}) 

class(sweden.c)
## [1] "sf"         "data.frame"
### water data ----
input_vector
## [1] "input_data/vector/animal_data/data_cleaned_PV.rds"                                      
## [2] "input_data/vector/animal_data/data_ready_rsf.rds"                                       
## [3] "input_data/vector/animal_data/reindeer_data.rds"                                        
## [4] "input_data/vector/country/swe_no.shp"                                                   
## [5] "input_data/vector/country/swe.shp"                                                      
## [6] "input_data/vector/HR/spatial_files_availability_PV_winter_2021_shorter_modified_new.shp"
## [7] "input_data/vector/ne_10m_lakes/ne_10m_lakes.shp"
# import lakes
lakes <-  input_vector %>%
  grep(pattern = "lakes", value = T) %>%
  sf::st_read() %>%
  st_as_sf() 
## Reading layer `ne_10m_lakes' from data source 
##   `/Users/hirn0001/Library/CloudStorage/OneDrive-Sverigeslantbruksuniversitet/Undervisning/Teaching and courses/GS_VMAS Data handling and illustrations/2025/create-maps/input_data/vector/ne_10m_lakes/ne_10m_lakes.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1355 features and 41 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -165.9656 ymin: -50.66967 xmax: 177.1544 ymax: 81.95521
## Geodetic CRS:  WGS 84
class(lakes)
## [1] "sf"         "data.frame"
st_crs(lakes)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
# Transform or convert coordinates of simple feature (sf)
lakes <- st_transform(lakes, crs = 3006) 
suppressWarnings({ 
  
# crop an sf object to a specific rectangle (here based on extent of "sweden.c"-shp)
# sf_use_s2(FALSE) # use if errors occur below
lakes <- st_crop(lakes, sweden.c) 

class(lakes)
plot(lakes)
}) 

Create points of interest (cities)

Create xy points with labels. Here creating sf-object and setting the crs.

# Create points of interest (coordinates from QGIS or google maps for demonstration)
places <- data.frame(ID= c("Arjeplog", "Arvidsjaur", "Sorsele"),
                     y = c(7328407.945, 7280856.153, 7270326.684),
                     x = c(630040.101, 693262.414, 617477.638)) %>%
  
  st_as_sf(coords = c("x", "y"),  # create sf object 
           
           crs=st_crs(sweden.c))  # set same crs as for sweden.c     
class(places)
## [1] "sf"         "data.frame"
 st_crs(places) # check projection 
## Coordinate Reference System:
##   User input: SWEREF99 TM 
##   wkt:
## PROJCRS["SWEREF99 TM",
##     BASEGEOGCRS["SWEREF99",
##         DATUM["SWEREF99",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4619]],
##     CONVERSION["SWEREF99 TM",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",15,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Topographic mapping (medium and small scale)."],
##         AREA["Sweden - onshore and offshore."],
##         BBOX[54.96,10.03,69.07,24.17]],
##     ID["EPSG",3006]]

Make plot using ggplot

library(ggspatial)
plot(sweden.c)
## Warning: plotting the first 9 out of 168 attributes; use max.plot = 168 to plot
## all

plot(lakes)
## Warning: plotting the first 9 out of 41 attributes; use max.plot = 41 to plot
## all

plot(places)

ggplot() +
  
    # plot sweden and lakes 
    geom_sf(aes(geometry = sweden.c$geometry)) + # note: geom_sf() for vector data 
  
    geom_sf(data = sweden.c,                     # plot sweden border vector
          aes(geometry = sweden.c$geometry),
          fill= "antiquewhite")+

      geom_sf(data = lakes,                      # plot lakes vector
              fill="lightblue") +

  # layer_spatial(data=elevation) +
  
  # add cities 
  geom_sf(data = places) +
  
  # add city names 
  geom_sf_text(data = places, aes(label = ID),
               vjust=-0.5)+
  
  # Add text on map
  annotate(geom = "text", 
             x = 484609, y = 7404230, 
             label = "Norway", 
             fontface = "italic", 
             color = "darkgray", 
           size = 6) +
  
    annotate(geom = "text", 
             x = 784609, y = 7404230, 
             label = "Sweden", 
             fontface = "italic", 
             color = "darkgray", 
           size = 6) +
  
  # Add north arrow
  annotation_north_arrow(location = "bl", 
                       which_north = "true", 
                         height = unit(1, "cm"),
  width = unit(1, "cm"),
                       pad_x = unit(0.5, "cm"), # horizontal align
                       pad_y = unit(8, "cm"), # vertical align
                       style = north_arrow_fancy_orienteering) +
  #add annotation scale
  annotation_scale(location = "bl", 
                     width_hint = 0.5) +
  
  
  # keep new projection for ggplot
  coord_sf(crs = st_crs(3006), 
  datum = sf::st_crs(3006),             
  expand = FALSE) +
  
  # add grid lines 
  theme(panel.grid.major = element_line(color = gray(.5), 
                                        linetype = "dashed", 
                                        linewidth  = 0.1), 
        panel.background = element_rect(fill = "lightblue")) 

### examplefor removing or adding titles 
# theme(axis.title.x = element_blank(),    # remove 
#       axis.title.y = element_blank()) + 
  # xlab("Longitude") + ylab("Latitude")   # add

ggsave("output_data/study_area.png", 
       width = 8, 
       height = 6, 
       dpi = 400)

Optional: Create a polygon and add to map

Example 2 interactive maps using mapview

# some gps data from reindeer 
reindeer_data <- readRDS("input_data/vector/animal_data/reindeer_data.rds") %>%
   dplyr::select(Collar_ID,  id, year, x_, y_, t_, group) 

str(reindeer_data)
## 'data.frame':    5951 obs. of  7 variables:
##  $ Collar_ID: chr  "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
##  $ id       : chr  "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
##  $ year     : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ x_       : num  733706 733653 733836 733828 733683 ...
##  $ y_       : num  7314782 7314737 7314807 7314851 7314738 ...
##  $ t_       : chr  "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
##  $ group    : chr  "Fed" "Fed" "Fed" "Fed" ...
### basic usage ===========================================================

### most basic call
mapview()
### with some data
# mapview(reindeer_data) 
str(reindeer_data)
## 'data.frame':    5951 obs. of  7 variables:
##  $ Collar_ID: chr  "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
##  $ id       : chr  "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
##  $ year     : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ x_       : num  733706 733653 733836 733828 733683 ...
##  $ y_       : num  7314782 7314737 7314807 7314851 7314738 ...
##  $ t_       : chr  "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
##  $ group    : chr  "Fed" "Fed" "Fed" "Fed" ...
# convert to sf 
reindeer_data_sf <- reindeer_data %>%
  st_as_sf(coords = c("x_", "y_"), #data frame to sf 
           crs=3006)              # assign to known crs 



plot(reindeer_data_sf)

mapview(reindeer_data_sf) 
# mapview structure
m <- mapview(reindeer_data_sf) 
m